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We investigate deterministic diffusion in periodic billiard models, in terms of the convergence of rescaled dis- 
tributions to the limiting normal distribution required by the central limit theorem; this is stronger than the usual 
requirement that the mean square displacement grow asymptotically linearly in time. The main model studied 
is a chaotic Lorentz gas where the central limit theorem has been rigorously proved. We study one-dimensional 
position and displacement densities describing the time evolution of statistical ensembles in a channel geometry, 
using a more refined method than histograms. We find a pronounced oscillatory fine structure, and show that 
this has its origin in the geometry of the billiard domain. This fine structure prevents the rescaled densities from 
converging pointwise to gaussian densities; however, demodulating them by the fine structure gives new densi- 
ties which seem to converge uniformly. We give an analytical estimate of the rate of convergence of the original 
distributions to the limiting normal distribution, based on the analysis of the fine structure, which agrees well 
with simulation results. We show that using a Maxwellian (gaussian) distribution of velocities in place of unit 
speed velocities does not affect the growth of the mean square displacement, but changes the limiting shape of 
the distributions to a non-gaussian one. Using the same methods, we give numerical evidence that a non-chaotic 
polygonal channel model also obeys the central limit theorem, but with a slower convergence rate. 

PACS numbers: 05.45.Pq, 02.50.-r, 02.70.Rr, 05.40.Jc 



I. INTRODUCTION 

Diffusion, the process by which concentration gradients are 
smoothed out, is one of the most fundamental mechanisms in 
physical systems out of equilibrium. Understanding the mi- 
croscopic processes which lead to diffusion on a macroscopic 
scale is one of the goals of statistical mechanics [1]. Since 
Einstein's seminal work on Brownian motion [2], diffusion 
has been modeled by random processes. However, we ex- 
pect the microscopic dynamics to be described by determinis- 
tic equations of motion. 

Recently it has been realized that many simple determin- 
istic dynamical systems are diffusive in some sense; we call 
this deterministic diffusion. Such systems can be regarded as 
toy models to understand transport processes in more realistic 
systems [1]. Examples include classes of uniformly hyper- 
bolic one-dimensional (ID) maps (see e.g. [3] and references 
therein) and multibaker models [4]. Often rigorous results are 
not available, but numerical results and analytical arguments 
indicate that diffusion occurs, for example in hamiltonian sys- 
tems such as the standard map [5]. 

Billiard models, where non-interacting point particles in 
free motion undergo elastic collisions with an array of fixed 
scatterers, have been particularly studied, since they are re- 
lated to hard sphere fluids, while being amenable to rigorous 
analysis [4, 6, 7]. They can also be regarded as the simplest 
physical systems in which diffusion, understood as the large- 
scale transport of mass through the system, can occur [8]. In 
this paper we study deterministic diffusion in two 2D billiard 
models: a periodic Lorentz gas, where the scatterers are dis- 
joint disks, and a polygonal billiard channel. 

A definition often used in the physical literature is that a 
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system is diffusive if the mean square displacement grows 
proportionally to time t , asymptotically as t — > oo. However, 
there are stronger properties which are also characteristic of 
diffusion, which a given system may or may not possess: (i) 
a central limit theorem may be satisfied, i.e. rescaled distribu- 
tions converge to a normal distribution as t — ► °°; and (ii) the 
rescaled dynamics may 'look like' Brownian motion. 

Two-dimensional (2D) periodic Lorentz gases were proved 
in [6, 7] to be diffusive in these stronger senses if they satisfy 
a geometrical finite horizon condition (Sec. II A). We use a 
square lattice with an additional scatterer in each cell to sat- 
isfy this condition, a geometry previously studied in [9, 10]. 
This model is of interest since, unlike in the commonly stud- 
ied triangular lattice case (see e.g. [4, 11, 12]), we can vary in- 
dependently two physically relevant quantities: the available 
volume in a unit cell, and the size of its exits; this is possible 
due to the two-dimensional parameter space [13, 14]. 

The main focus of this paper is to investigate the fine struc- 
ture occurring in the position and displacement distributions 
at finite time t , and the relation with the convergence to a lim- 
iting normal distribution as t — > °° proved in [6, 7]. Those 
papers show in what sense we can smooth away the fine struc- 
ture to obtain convergence. However, from a physical point 
of view it is important to understand how this convergence 
occurs; our analysis provides this. 

This analysis makes explicit the obstruction that prevents 
a stronger form of convergence, showing how density func- 
tions fail to converge pointwise to gaussian densities; it also 
allows us to conjecture a more refined result which takes the 
fine structure into account. 

Furthermore, this line of argument suggests how conver- 
gence may occur in other models where few rigorous re- 
sults are available. As an example, we analyze a recently- 
introduced polygonal billiard channel model, showing that the 
same techniques are still applicable. 



Plan of paper 

In Sec. II we present the periodic Lorentz gas model for 
which we obtain most of our results. Section III discusses the 
definition of diffusion in the context of deterministic dynam- 
ical systems. In Sec. IV we study numerically the fine struc- 
ture of distributions in the Lorentz gas, finding good agree- 
ment with an analytical calculation in terms of the geometry of 
the billiard domain, and showing that when this fine structure 
is removed, the demodulated densities are close to gaussian. 
This we apply in Sec. V to investigate the central limit theo- 
rem and the rate of convergence to the limiting normal distri- 
bution, obtaining a simple estimate of this rate which agrees 
well with numerical results. In Sec. VI we study the effect 
of imposing a Maxwellian (gaussian) velocity distribution in 
place of a unit speed distribution, showing that this leads to 
non-gaussian limiting distributions. Section VII extends these 
ideas to a polygonal billiard channel, where few rigorous re- 
sults are available. We finish with conclusions in Sec. VIII. 
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FIG. 1 : (a) Part of the infinite system, constructed from two square 
lattices of disks shown in different shades of gray; dashed lines indi- 
cate several unit cells and an elastic collision is shown, (b) A single 
unit cell, defining the geometrical parameters. The billiard domain is 
the area Q exterior to the disks. 



n. TWO-DIMENSIONAL PERIODIC LORENTZ GAS 

We consider periodic billiard models, where the dynamics 
can be studied on the torus. The region Q exterior to the scat- 
terers is called the billiard domain; we denote its area by \Q\. 
Since the particles are non-interacting, it is usual to set all ve- 
locities to 1 by a geometrical rescaling, although in Sec. VI 
we discuss the effect of a gaussian velocity distribution. 

We focus on a periodic Lorentz gas, where the scatterers are 
non-overlapping disks. Their strictly convex boundaries make 
this a scattering billiard [6], and hence a chaotic system, in 
the sense that it has a positive Lyapunov exponent [4, 15] and 
positive Kolmogorov-Sinai entropy [4]. 



A. Periodic Lorentz gas model 

The model we study, previously considered in [9, 10], con- 
sists of two square lattices of disks; they have the same lattice 
spacing r, and radii a and b, respectively, and are positioned 
such that there is a fe-disk at the center of each unit cell of 
the a-lattice: see Fig. 1. In analytical calculations we take the 
length scale as r — 1, as in [9, 10], whereas in numerical sim- 
ulations we fix a = 1 and scale r and b appropriately, as in 
[12]. 

Finite horizon condition Periodic Lorentz gases were 
shown in [6, 7] to be diffusive (Sec. Ill), provided they sat- 
isfy the finite horizon condition: there is an upper bound on 
the free path length between collisions. If this is not the case, 
so that a particle can travel infinitely far without colliding 
(the billiard has an infinite horizon), then corridors exist [16], 
which allow for fast propagating trajectories, leading to super- 
diffusive behavior, as was recently rigorously proved [17]. 

We restrict attention to parameter values within the finite 
horizon regime by choosing b to block all corridors [10, 13]. 



B. Statistical properties 

Statistical properties of deterministic dynamical systems 
arise from an ensemble of initial conditions (xq,Vo) model- 
ing the imprecision of physical measurements. We always 
take a uniform distribution with respect to Liouville measure 
in one unit cell: the positions xo are uniform with respect to 
Lebesgue measure in the billiard domain Q, and the velocities 
Vo are uniform in the unit circle S , i.e. with angles between 
and 2k, and unit speeds. 

We evolve (xn,vn) for a time t under the billiard flow <t>' 
in phase space to (x(f),v(f)). Note that Liouville measure 
on the torus is invariant under this flow [15]. In numerical 

experiments, we take a large sample (x$\ v$)%=\ °f s i ze N °f 
initial conditions chosen uniformly with respect to Liouville 
measure using a random number generator. These evolve after 
time t to {*('\t)M'\t))f =l ; the distribution of this ensemble 
then gives an approximation to that of (x(r), v(t )). 

We denote averages over the initial conditions, or equiva- 
lently expectations with respect to the distribution of (xq, v ), 
by (•). Approximations of such averages can be evaluated us- 
ing a simple Monte Carlo method [18] as 




(2.1) 



The infinite sample size limit, although unobtainable in prac- 
tice, reflects the expectation that larger N will give a better 
approximation. Averages at time t can be evaluated by using 
a function / involving <t>' . 

C. Channel geometry 

Diffusion occurs in the extended system obtained by un- 
folding the torus to a 2D infinite lattice: see [6, 7] and Sec. III. 
The diffusion process is then described by a second order dif- 
fusion tensor having 4 components Dy with respect to a given 
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FIG. 2: ID channel obtained by unfolding torus in x-direction. 
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FIG. 3: (a) Lorentz channel studied in [19, 20] with hard upper and 
lower boundaries; dotted lines indicate unit cells, (b) Fully unfolded 
triangular Lorentz gas. Dotted lines indicate unit cells forming a 
channel with periodic upper and lower boundaries. 



A. Diffusion as a stochastic process 

Diffusion is described classically by the diffusion equation 
^M=DV^( f ,x), (3.1) 

where p is the density of the diffusing substance. Following 
Einstein and Wiener (see e.g. [2]), we can model diffusion as 
a stochastic process B t , determined by the probability density 
p(x,t) of a particle being at position x at time t given that it 
started at x = at time t = 0. 

Imposing conditions on the process determined from phys- 
ical requirements gives a diffusion process, where p(x,t) sat- 
isfies the equation 



dp_ d_ 

dt dx{ 



= o, 



(3.2) 



known as Kolmogorov's forward equation or the Fokker- 
Planck equation [2]. The drift vector A(x,f) and the diffu- 
sion tensor B(x,f) give the mean and variance, respectively, 
of infinitesimal displacements at position x and time t [2]. 

If the system is sufficiently symmetric that the drift is zero 
and the diffusion tensor is a multiple of the identity tensor, 
then the process is Brownian motion, and (3.2) reduces to the 
diffusion equation (3.1). A general diffusion process, how- 
ever, can be inhomogeneous in both space and time. 



orthonormal basis, given by 

D ij = lim^-(Ax i Ax j ) t . (2.2) 
2f 

The square symmetry of our model reduces the diffusion 
tensor to a constant multiple D of the identity tensor; we can 
evaluate this diffusion coefficient by restricting attention to 
the dynamics in a 1 -dimensional channel extended only in the 
x-direction; see Fig. 2. Correspondingly, we restrict attention 
to ID marginal distributions. 

A channel geometry, with hard horizontal boundaries, cor- 
responding to the triangular Lorentz gas was studied in [19, 
20] (Fig. 3(a)). This is equivalent to a channel with twice the 
original height and periodic boundaries, shown in Fig. 3(b) 
as part of the whole triangular lattice obtained by unfolding 
completely in the vertical direction. We can view this lattice 
as consisting of rectangular unit cells (Fig. 3(b)) which are 
stretched versions of the square unit cell considered above, 
with the extra condition a = b. The results in the remainder of 
this paper then extend to this case with minor changes. 



III. DETERMINISTIC DIFFUSION 

In this section we briefly recall how to make precise the fact 
that the behavior of certain deterministic dynamical systems 
'looks like' that of the diffusion equation. 



B. Diffusion in dynamical systems via limit theorems 

Diffusion in billiards concerns the statistical behavior of the 
particle positions. We can write the first component x t of the 
position x, at time t as 

x,= [ Vl (s)ds+x = f fo<t> s (-)ds+x Q , (3.3) 
Jo Jo 

where / = vi, the first velocity component. This expresses x t 
solely in terms of functions defined on the torus. In fact, (3.3) 
shows that the displacement Ax t := x t — xq is in some sense a 
more natural observable than the position x t in this context. 

We thus wish to study the distribution of accumulation 
functions of the form S t {-) :— Jq f o <P S (■) ds, in particular in 
the limit as t — > °° [21]. We remark that other observables / 
are relevant for different transport processes [8]. 

We denote by <t>' : — > the flow of a dynamical system 
with time t £ K. Given a probability measure fi describing the 
distribution of initial conditions, we can find the probability of 
being in certain regions of the phase space ^ at given times, 
so that we have a stochastic process. If the measure jtt is in- 
variant, so that jU(4>~'(A)) = ju(A) for all times t and all nice 
sets A, then the stochastic process is stationary [21]. 

The integral in the definition of S is then a continuous- 
time version of a Birkhoff sum Y.1Zo f ° over tne station- 
ary stochastic process given by <t>, so that we may be able to 
apply limit theorems from the theory of stationary stochastic 
processes [21]. For the case of the periodic Lorentz gas with 
finite horizon, it was proved in [6, 7] that the following limit 
theorems hold. 
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a. Asymptotic linearity of mean square displacement 
The limit 

2D := Mm -(Ax 2 ), (3.4) 

exists, so that the mean square displacement (Ax 2 ), := 
([Ax(f)] 2 ) (the variance of the displacement distribution) 
grows asymptotically linearly in time: 

(Ac 2 ) ( ~ 2Dt as t -> °°, (3.5) 

where D is the diffusion coefficient. In d > 2 dimensions, 
setting Axj(t) := x,(f) — x,(0), we have 

(AxiAxj),~2Dijt, (3.6) 

where the are components of a symmetric diffusion tensor. 

b. Central limit theorem: convergence to normal distri- 
bution Scale the displacement distribution by y/i, so that the 
variance of the rescaled distribution is bounded. Then this dis- 
tribution converges weakly, or in distribution, to a normally 
distributed random variable z [21, 22]: 

x(f)-x(0) 9 

' — ► z, as t — > oo. (3.7) 

v f 

In the 1 -dimensional case, this means that 

HmP M <H ) 1 f e-^ 2 dv , (3.8) 

where P(-) denotes probability with respect to the distribution 
of the initial conditions, and a 2 is the variance of the limiting 
normal distribution. In d > 2 dimensions, this is replaced by 
similar statements about probabilities of <f-dimensional sets. 
This is the central limit theorem for the random variable Ax. 
From (a) we know that in ID, the variance of the limiting nor- 
mal distribution is a 2 = 2D; \nd>2 dimensions, the covari- 
ance matrix of z is given by the matrix (2D,-;) [6, 23]. 

c. Functional central limit theorem: convergence of path 
distribution to Brownian motion We rescale the path x t by 
the scale from (b), defining % by [16] 

M s):= X -^m, ,e[0,l]. (3.9) 

vt 

The distribution of these rescaled paths then converges in dis- 
tribution to Brownian motion: 

x f ^B asf^oo, (3.10) 

where the Brownian motion B has co variance matrix as in (b). 
This is known as a functional central limit theorem, or weak 
invariance principle [21]. 

A sufficient condition for this is that the following two 
properties hold [24]. (i) The multi-dimensional central limit 
theorem, a generalization of (b), is satisfied. This says that 
the finite-dimensional distributions of the process x t converge 
to those of Brownian motion, so that for any n, any times 



s\ < ■ ■ ■ < s n , and any reasonable sets D\,... ,D„ in R d , we 
have 

V(x t {s 1 )eD l ,...,x t (s n )eD n ) 

P (B(ji) eDi,... ,B(sn) G A,) • (3.11) 

The right-hand side can be expressed as a multi-dimensional 
integral overgaussians: see e.g. [14, 23]. (ii) The convergence 
is tight, which prevents mass escaping to infinity: see [24] for 
the definition. 



C. Discussion of definitions of diffusion 

Property (c) is the strongest sense in which a dynamical 
system can show deterministic diffusion, making precise how 
a rescaled dynamical system can look like Brownian motion. 
However, few physically relevant systems have been proved to 
satisfy (c): interest in the periodic Lorentz gas comes largely 
from the fact that it is one; another is the triple linkage [25]. 

The multi-dimensional central limit theorem part of (c) was 
studied in [23], where both Lorentz gases and wind-tree mod- 
els were found to obey it, tested for certain sets D, and certain 
values of n. However, as stated in [23], (c) is difficult to in- 
vestigate numerically, and the results in that paper seem to be 
the best that we can expect. 

Property (b), the central limit theorem, has been shown for 
large classes of observables / in many dynamical systems (see 
[21] and references therein), but again they are often not phys- 
ical. Property (b) was used in [22] as the definition of a dif- 
fusive system, but does not seem to have been applied in the 
physical literature; it is the approach taken in this paper. 

Many papers in the physical literature define a system to 
be diffusive if only property (a) is verified (numerically), e.g. 
[12, 26, 27]. Many types of system are diffusive in this sense, 
including ID maps [3], random Lorentz gases [27] and Ehren- 
fest wind-tree models, both periodic [26] and random [27]. 

It is possible for the weaker properties to hold when the 
stronger ones do not. For example, in [28] a disordered lattice- 
gas wind-tree model was reported to have an asymptotically 
linear mean square displacement, but a non-gaussian distri- 
bution function, i.e. (a) but not (b). However, disorder can 
lead to trapping effects which cannot occur in periodic sys- 
tems [26], and we are not aware of a periodic (and hence or- 
dered) billiard-type model with unit-speed velocity distribu- 
tion which shows (a) but not (b), although in Sec. VI we show 
that this can occur with a Maxwellian velocity distribution. 

IV. FINE STRUCTURE OF POSITION AND 
DISPLACEMENT DISTRIBUTIONS 

We now focus on the diffusive properties of the periodic 
Lorentz gas model introduced in Sec. II. In this section, we 
describe the fine structure of position and displacement distri- 
butions. The displacement distribution occurs naturally in the 
central limit theorem (Sec. Ill B) and in Green-Kubo relations 
[1,4], whereas the position distribution is more natural if we 
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are unable to track the paths of individual particles. It is possi- 
ble to show that the asymptotic properties of the position and 
displacement distributions are the same, in the sense that one 
has an asymptotically linear growth if and only if the other 
does, and similarly for the central limit theorem [14]. It is 
hence equivalent to consider diffusive properties by studying 
either distribution. 



A. Position and displacement distributions 

Figure 4 shows scatterplots representing 2D position and 
displacement distributions for a representative choice of geo- 
metrical parameters. Each dot represents one initial condition 
started in the central unit cell and evolved for time t = 50; 
N = 5 x 10 4 samples are shown. Both distributions show de- 
cay away from a maximum in the central cell, an overall cir- 
cular shape, and the occurrence of a periodic fine structure. 

These figures are projections to the billiard domain Q of 
the density in the phase space Q x S 1 . Since the dynamics on 
the torus is mixing [15], the phase space density converges 
weakly [29] to a uniform density on phase space correspond- 
ing to the invariant Liouville measure. Physically, the phase 
space density develops a complicated layer structure in the 
stable direction of the dynamics: see e.g. [1]. Projecting cor- 
responds to integrating over the velocities; we expect this to 
eliminate this complicated structure and result in some degree 
of smoothness of the projected densities. However, we are not 
aware of any rigorous results in this direction, even for rel- 
atively well-understood systems such as the Arnold cat map 
[1]. 

These 2D distributions are difficult to work with, and we 
instead restrict attention to one-dimensional marginal distri- 
butions, i.e. projections onto the x-axis, which will also have 
some degree of smoothness. We denote the ID position den- 
sity at time t and position x € K by f t (x) and the displacement 
density for displacement x by gt(x). We let their respective 
(cumulative) distribution functions be F t (x) and G t (x), respec- 
tively, so that 

F t (x) := P(x t < x) = f f t (s)ds, (4.1) 

./ — oo 

and similarly for G t . (When necessary, we will instead denote 
displacements by £, .) The densities show the structure of the 
distributions more clearly, while the distribution functions are 
more directly related to analytical considerations. 



B. Numerical estimation of distribution functions and densities 

We wish to estimate numerically the above denstities 
and distribution functions at time t from the N data points 
x\ l \... ,Xf N \ The most widely used method in the physics 
community for estimating density functions from numerical 
data is the histogram; see e.g. [26]. However, histograms are 
not always appropriate, due to their non-smoothness and de- 
pendence on bin width and position of bin origin [30]. In [26], 



for example, the choice of a coarse bin width obscured the fine 
structure of the distributions that we describe in Sec. VII. 

We have chosen the following alternative method, which 
seems to work well in our situation, since it is able to deal with 
strongly peaked densities more easily, although we do not 
have any rigorous results to justify it. We have also checked 
that histograms and kernel density estimates (a generalization 
of the histogram [30]) give similar results, provided sufficient 
care is taken with bin widths. 

We first calculate the empirical cumulative distribution 

function [30, 31], defined by F f emp (x) := #{; : x^ < x] for 
the position distribution, and analogously for the displacement 
distribution. The estimator F t emp is the optimal one for the 
distribution function F, given the data, in the sense that there 
are no other unbiased estimators with smaller variance [31, p. 
34]. We find that the distribution functions in our models are 
smooth on a scale larger that that of individual data points, 
where statistical noise dominates. (Here we use 'smooth' in 
a visual, nontechnical sense; this corresponds to some degree 
of differentiability). We verify that adding more data does not 
qualitatively change this larger-scale structure: with N =\Q 1 
samples we seem to capture the fine structure. 

We now wish to construct the density function f t = dF t jdx. 
Since the direct numerical derivative of F t emp is useless due to 
statistical noise, our procedure is to fit an (interpolating) cu- 
bic spline to an evenly-spread sample of points from F t emp , 
and differentiate the cubic spline to obtain the density func- 
tion at as many points as required [14]. Sampling evenly 
from F, cmp automatically uses more samples where the data 
are more highly concentrated, i.e. where the density is larger. 

We must confirm (visually or in a suitable norm) that our 
spline approximation reproduces the fine structure of the dis- 
tribution function sufficiently well, whilst ignoring the vari- 
ation due to noise on a very small scale. As with any den- 
sity estimation method, we have thus made an assumption of 
smoothness [30]. The analysis of the fine structure in Sec. IV 
justifies this to some extent. 



C. Time evolution of ID distributions 

Figure 5 shows the time evolution of ID displacement dis- 
tribution functions and densities for certain geometrical pa- 
rameters, chosen to emphasize the oscillatory structure. Other 
parameters within the finite horizon regime give qualitatively 
similar behavior. 

The distribution functions are smooth, but have a step-like 
structure. Differentiating the spline approximations to these 
distribution functions gives densities which have an under- 
lying gaussian-like shape, modulated by a pronounced fine 
structure which persists at all times (Fig. 5(b)). This fine 
structure is just noticeable in Figs. 4 and 5 of [26], but other- 
wise does not seem to have been reported previously, although 
in the context of iterated ID maps a fine structure was found, 
the origin of which is pruning effects: see e.g. Fig. 3.1 of [32]. 
We will show that in billiards this fine structure can be under- 
stood by considering the geometry of the billiard domain. 
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FIG. 4: (a) 2D position distribution; (b) 2D displacement distribution, r = 2.5; b = 0.4; t = 50; JV = 5x 10 initial conditions. 




FIG. 5: (Color online) (a) Time evolution of displacement distribution functions, (b) Time evolution of displacement densities, calculated by 
numerically differentiating a cubic spline approximation to distribution functions, r = 2.1; b = 0.2. The inset in (a) shows the definition of the 
set H(x) required later. 



D. Fine structure of position density 

Since Liouville measure on the torus is invariant, if the ini- 
tial distribution is uniform with respect to Liouville measure, 
then the distribution at any time t is still uniform. Integrating 
over the velocities, the position distribution at time t is hence 
always uniform with respect to Lebesgue measure in the bil- 
liard domain Q, which we normalize such that the measure of 
Q is 1. Denote the two-dimensional position density on the 
torus at (x,y) £ [0, l) 2 by p toms (x,y). Then 



121 



)(y). 



(4.2) 



Here, H(x) := {y : (x,y) £ Q} is the set of allowed y values 
for particles with horizontal coordinate x (Fig. 5(a) inset), and 
1b is the indicator function of the (one- or two-dimensional) 
set B, given by 



U{b) 



1, \fb£B 
0, otherwise. 



(4.3) 



Thus for fixed x, p toms (x,y) is independent of y within the 
available space H(x). 

Now unfold the dynamics onto a 1 -dimensional channel in 
the x-direction, as in Fig. 2, and consider the torus as the 
distinguished unit cell at the origin. Fix a vertical line with 
horizontal coordinate x in this cell, and consider its periodic 
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translates x + n along the channel, where n G Z. Denoting the 
density there by p r channel (x + n,y), we have that for all t and 
for all x and y, 

^p,"(x + n,y) = p tOTUS (x,y). (4.4) 

We expect that after a sufficiently long time, the distribu- 
tion within cell n will look like the distribution on the torus, 
modulated by a slowly-varying function of x. In particular, 
we expect that the 2D position density will become asymptot- 
ically uniform in y within H(x) at long times. We have not 
been able to prove this, but we have checked by constructing 
2D kernel density estimates [30] that it seems to be correct. 
A 'sufficiently long' time would be one which is much longer 
than the time required for the diffusion process to cross one 
unit cell. 

Thus we have approximately 

p,"(x,y) ~ p toms (x,y) p t (x) = ft (,)-Ll flW (y), (4.5) 

where pt{x) is the shape of the two-dimensional density dis- 
tribution as a function of x G R; we expect this to be a slowly- 
varying function. We use '~' to denote that this relationship 
holds in the long-time limit, for values of x which do not lie in 
the tails of the distribution. Although this breaks down in the 
tails, the density is in any case small there. 

The ID marginal density that we measure will then be given 
approximately by 

ft{x)=f pf"(x,y)dy^p t (x)h(x), (4.6) 

Jy=0 

where h(x) := |i/(jt)| / \Q\ is the normalized height (Lebesgue 
measure) of the set H(x) at position x (see the inset of 
Fig. 5(a)). Note that H(x) is not necessarily a connected set. 

Thus the measured density f t (x) is given by the shape p t (x) 
of the 2D density, modulated by fine-scale oscillations due to 
the geometry of the lattice and described by h(x), which we 
call the fine structure function. 

The above argument motivates the (re-)definition of p t (x) 
so that that f t (x) = h(x)p t {x), now with strict equality and for 
all times. We can then view p t (x) as the density with respect 
to a new underlying measure hX, where X is 1 -dimensional 
Lebesgue measure; this measure takes into account the avail- 
able space, and is hence more natural in this problem. We 
expect that p t will now describe the large-scale shape of the 
density, at least for long times and x comparatively small. 

Figure 6 shows the original and demodulated densities f t 
and p t for a representative choice of geometrical parameters. 
The fine structure in f t is very pronounced, but is eliminated 
nearly completely when demodulated by dividing by the fine 
structure h, leaving a demodulated density p, which is close 
to the gaussian density with variance 2Dt (also shown). 

We estimated the diffusion coefficient D as follows. For 
r = 2.3 and b = 0.5, using N = 10 7 particles evolved to 
t = 1000, the best fit line for log(Ax 2 ), against logf in the 
region t S [500, 1000] gives (Ax 2 ) ~ ? i oooo3^ wn j cn we regard 
as confirmation of asymptotic linear growth. Following [12], 





x 



FIG. 6: (Color online) Position density /, exhibiting a pronounced 
fine structure, together with the demodulated slowly-varying func- 
tion pi and a gaussian with variance 2Dt. The inset shows one pe- 
riod of the demodulating fine structure function h. r = 2.3; b = 0.5; 
t = 50. 

we use the slope of log(Ax 2 ), against t in that region as an 
estimate of 2D, giving D = 0.1494 ±0.0002; see [13, 14] for 
the error analysis. 

(Throughout the paper, we denote by y a 2 the gaussian den- 
sity with mean and variance rj 2 , and by N a 2 the correspond- 
ing normal distribution function.) 

Note that although the density has non-smooth points, 
which affects the smoothness assumption in our density es- 
timation procedure described in Sec. IV B, in practice these 
points are still handled reasonably well. If necessary, we could 
treat these points more carefully, by suitable choices of parti- 
tion points in that method. 



E. Fine structure of displacement density 

We can treat the displacement density similarly, as follows. 
Let rj t (x,y) be the 2D displacement density function at time t, 
so that 

J j\t(x,y)dxdy = F(Ax t <x,Ay, <y), x,yeR. 

(4.7) 

(Recall that Ax t := x t —xq.) We define the projected versions 

^channel an( j ^torus afj follows: 

r/ «(x,y):=£T},(x,y + n), xeR,ye[0,l), (4.8) 
nr s (x,y) := £ 7 7 r Channel (^ + ",3') I x,y E [0, 1). (4.9) 

Again we view the torus as the unit cell at the origin where all 
initial conditions are placed. Note that projecting the displace- 
ment distribution on R 2 to the channel or torus gives the same 
result as first projecting and then obtaining the displacement 
distribution in the reduced geometry. Hence the designations 
as being associated with the channel or torus are appropriate. 
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Unlike p torus in the previous section, rj r toms is not indepen- 
dent of t : for example, for small enough t , all displacements 
increase with time. However, we show that Tj f torus rapidly ap- 
proaches a distribution which is stationary in time. 

Consider a small ball of initial conditions of positive Liou- 
ville measure around a point (x, v). Since the system is mixing 
on the torus, the position distribution at time t corresponding 
to those initial conditions converges as t — > °° to a distribu- 
tion which is uniform with respect to Lebesgue measure in the 
billiard domain Q. The corresponding limiting displacement 
distribution is hence obtained by averaging the displacement 
of x from all points on the torus. 

Extending this to an initial distribution which is uniform 
with respect to Liouville measure over the whole phase space, 
we see that the limiting displacement distribution is given by 
averaging displacements of two points in Q, with both points 
distributed uniformly with respect to Lebesgue measure on Q. 
This limiting distribution we denote by 77 torus (x,y), with no t 
subscript. 

As in the previous section, we expect the y-dependence of 



channel 



(x + n, •) to be the same, for large enough t, as that of 



?7 toms (x, ^ f or x e [o, 1). However, 77 toms (x, •) is not indepen- 
dent of >', as can be seen from a projected version of Fig. 4(b) 
on the torus [14]. We thus set 



, channel 



(x,y)^7 toras (x,y)T7 f (x). 



(4.10) 



To obtain the ID marginal density g t {x), we integrate with 
respect to y: 



8t(x) 



f 

Jy= 



, channel 



(x,y)dy~ <j>(x)i) t {x), 



where 



0(x) := 



Jy=0 



5 (x,y)dy. 



(4.12) 



Again we now redefine fj so that gt(x) — (p(x)fj t (x), with the 
fine structure of g t (x) being described by and the large-scale 
variation by fj(x), which can be regarded as the density with 
respect to the new measure X taking account of the excluded 
volume. In the next section we evaluate (x) explicitly. 



This form leads us to expand in Fourier series: 

h(x)= Y d Hk)s 2Kikx ^h(^)+2 1 £kk)cos2nkx, (4.14) 

kez keN 

and similarly for 0, where the Fourier coefficients are defined 
by 

h(k):= J h(x)e- 2mkx dx = J h(x) cos(2xkx)dx. (4.15) 

The last equality follows from the evenness of h, and shows 
that h(k) = h(—k), from which the second equality in (4.14) 
follows. Fourier transforming (4.13) then gives 



§(k)=h(k)h{-k)=h{kf 



(4.16) 



Taking the origin in the center of the disk of radius b (see 
the inset of Fig. 5), the available space function h is given by 



h(x) 



1 

W\ 



\-2\/b 2 -x 2 -2 



(4.17) 



forx € [0, 1 /2), and is even and periodic with period 1. Here 
we adopt the convention that y/a = if a < to avoid writing 
indicator functions explicitly. The evaluation of the Fourier 
coefficients of h thus involves integrals of the form 



Jo 



coszf \/a 2 — t 2 dt = ^-Jiiza), 
2z 



(z^O) (4.18) 



where J\ is the first order Bessel function; this equality fol- 
lows from equation (9.1.20) of [33] after a change of vari- 
ables. 

(4-1 1) Hence the Fourier coefficients of h are h(0) = Jq h(x) = 1 



and, for integer k ^ 0, 

k{k) =~wm 



(-l) k aJ 1 (27ta\k\) +bJ 1 {2nb\k\) 



(4.19) 

Note that Jq <p (x) dx = (0) = h(Q) 2 = 1 , so that is correctly 
normalized as a density function on the torus. 

In Fig. 7 we plot partial sums m up to m terms of the 
Fourier series for analogous to (4.14). We can determine 
the degree of smoothness of 0, and hence presumably of gt, 
as follows. The asymptotic expansion of J\ (z) for large real z 
(equation (9.2.1) of [33]), 



F. Calculation of x-displacement density <j> (x) on torus 

Let (X\,Y\) and (X2,l2) be independent random variables, 
distributed uniformly with respect to Lebesgue measure in the 
billiard domain Q, and let AX := {X 2 -X{\ e [0, 1) be their 
x-displacement (where {•} again denotes the fractional part of 
its argument). Then AX is the sum of two independent ran- 
dom variables, so that its density is given by the following 
convolution, which correctly takes account of the periodicity 
of h and with period 1 : 



0(4)= / h(x)h(x + %)dx. 
Jo 



(4.13) 



/ 2 /37T 

/i(z) ~ \ — cos — — z 
V nz V 4 



= 0( z - l/2 



(4.20) 



shows thatfc(fc) = ^(^ 3/2 ) and hence 0(£) = <^(Ar 3 ). From 
the theory of Fourier series (see e.g. [34, Chap. 1]), we hence 
have that is at least C 1 (once continuously differentiable). 
Thus the convolution of h with itself is smoother than h is, as 
intuitively expected, despite the non-differentiable points of h. 

We have checked numerically the approach of 
f 77 ( torus (x, j y)d j y to 0(x), and it appears to be fast, al- 
though the rate is difficult to evaluate, since a large number 
of initial conditions are required for the numerically calcu- 
lated distribution function to approach closely the limiting 
distribution. 
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FIG. 7: (Color online) Partial sums <p,„ up to m terms of the Fourier 
series for <j), with r = 2.3 and ft = 0.5. 



G. Structure of displacement distribution 

In Fig. 8 we plot the numerically-obtained displacement 
density g t (x), the fine structure function (j) calculated above, 
and their ratio f\ t (x), for a certain choice of geometrical pa- 
rameters. Again the ratio is approximately gaussian, which 
confirms that the densities can be regarded as a gaussian shape 
modulated by the fine structure 0. 

However, if r is close to 2a, then fj, itself develops a type 
of fine structure: it is nearly constant over each unit cell. This 
is shown in Fig. 9 for two different times. We plot both g t and 
f) t , rescaled by \Jt and compared to a gaussian of variance 2D. 
(This scaling is discussed in Sec. V.) 

This step-like structure of fj t is related to the validity of the 
Machta-Zwanzig random walk approximation, which gives 
an estimate of the diffusion coefficient in regimes where the 
geometrical structure can be regarded as a series of traps with 
small exits [11-13, 35]. Having fj, constant across each cell 
indicates that the distribution of particles within the billiard 
domain in each cell is uniform, as is needed for the Machta- 
Zwanzig approximation to work. 

As r increases away from 2a, the exit size of the traps in- 
creases, and the Machta-Zwanzig argument ceases to give a 
good approximation [12, 13]. The distribution then ceases to 
be uniform in each cell: see Fig. 6. This may be related to the 
crossover to a Boltzmann regime described in [12]. 



V. CENTRAL LIMIT THEOREM AND RATE OF 
CONVERGENCE 

We now discuss the central limit theorem as t — > °° in terms 
of the fine structure described in the previous section. 



FIG. 8: (Color online) Displacement density g,, with demodulated 
fit compared to a gaussian of variance 2D. The inset in (a) shows the 
fine structure function (j) for these geometrical parameters, r — 2.1; 
fc = 0.2;f = 50. 




x 



FIG. 9: (Color online) Displacement density gt and demodulated fjt, 
both rescaled by \ft, at t — 200 and t = 1000, compared to a gaussian 
of variance 2D. The inset again shows the fine structure function (j). 
r = 2.01; b = 0.1. 



A. Central limit theorem: weak convergence to normal 
distribution 

The central limit theorem requires us to consider the densi- 
ties rescaled by \ft, so we define 

g,(x):=Vigt(xVt), (5.1) 

where the first factor of y/i normalizes the integral of g t to 1, 
giving a probability density. Figure 10 shows the densities of 
Fig. 5(a) rescaled in this way, compared to a gaussian density 
with mean and variance 2D. We see that the rescaled densi- 
ties oscillate within an envelope which remains approximately 
constant, but with an increasing frequency as t — > °o; they are 
oscillating around the limiting gaussian, but do not converge 
to it pointwise. See also Fig. 9. 

The increasingly rapid oscillations do, however, cancel out 
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A" 



FIG. 10: (Color online) Displacement densities as in Fig. 5(b) after 
rescaling by \ft, compared to a gaussian density with mean and 
variance 2D. r = 2.1; b = 0.2. 
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x 

FIG. 11: (Color online) Difference between rescaled distribution 
functions and limiting normal distribution with variance 2D. r = 2.1; 
b = 0.2. 



when we consider the rescaled distribution functions, given by 
the integral of the rescaled density functions: 

G,{x):= j*_ g t (s)6s = G t (xVt). (5.2) 

Figure 1 1 shows the difference between the rescaled distribu- 
tion functions and the limiting normal distribution with mean 
and variance 2D. We see that the rescaled distribution func- 
tions do converge to the limiting normal, in fact uniformly, as 
t — > oo; we thus have weak convergence. 

Although this is the strongest kind of convergence we can 
obtain for the densities g t with respect to Lebesgue measure, 
Fig. 9 provides evidence for the following conjecture: the 
rescaled densities f] t with respect to the new, modulated mea- 
sure converge uniformly to a gaussian density. This character- 
izes the asymptotic behavior more precisely than the standard 
central limit theorem. 
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FIG. 12: (Color online) Distance of rescaled distribution functions G, 
from limiting normal distribution A^d in log-log plot. The straight 
line is a fit to the large-time decay of the data for r = 2.05. 

B. Rate of convergence 

Since the G t converge uniformly to the limiting normal dis- 
tribution, we can consider the distance \\G t —N2o\\ oe , where 
we define the uniform norm by 

||F|U:=sup|F(x)|. (5.3) 

We denote by N a 2 the normal distribution function with vari- 
ance a 2 , given by 

N o2 (x):=-^= T e^^'ds, (5.4) 

Figure 12 shows a log-log plot of this distance against time, 
calculated numerically from the full distribution functions. 
We see that the convergence follows a power law 

\\&t-N2 D \\ Ba ~t- a , (5.5) 

with a fit to the data for r = 2.05 giving a slope a ~ 0.482. The 
same decay rate is obtained for a range of other geometrical 
parameters, although the quality of the data deteriorates for 
larger r, reflecting the fact that diffusion is faster, so that the 
distribution spreads further in the same time. Since we use the 
same number = 10 7 of initial conditions, there is a lower 
resolution near x — where, as shown in the next section, the 
maximum is obtained. 

In [36] it was proved rigorously that a > I ~ 0. 167 for any 
Holder continuous observable /. Here we have considered 
only the particular Holder observable v, but for this function 
we see that the rate of convergence is much faster than the 
lower bound proved in [36]. 

C. Analytical estimate of rate of convergence 

We now obtain a simple analytical estimate of the rate of 
convergence using the fine structure calculations in Sec. IV. 
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Since the displacement distribution is symmetric, we have 
G t (x = 0) = 1/2 for all t . The maximum deviation of G t from 
N2D occurs near to x = 0, where the density function is fur- 
thest from a gaussian, while the fine structure of the density g, 
means that G t is increasing there (Fig. 1 1). Due to the oscilla- 
tory nature of the fine structure, this maximum thus occurs at 
a distance of 1 /4 of the period of oscillation from x = 0. 

Since the displacement density has the form g t (x) = 
<j>(x)i) t (x), after rescaling we have 



g t (x) = <j)(xVt)rj(x), 



(5.6) 



where T] t [x):=\/i r] t (xy/t) is the rescaled slowly-varying part 
of g t , and the fine structure at time t is given by 

<t>{x\ft) = 1 + 2£ (j>(k) cos(2nkxVi). (5.7) 

ken 

The maximum deviation occurs at 1 /4 of the period of 
<j)(xy/t), i.e. at* = so that 

rl/VF 

\\Gt-NWn~ £ <P(k)cos(2nkxVt)dx (5.8) 

^_^(k-l)/2 



2nk 



(5.9) 



The correction due to the curvature of the underlying gaussian 
converges to as t — > », since this gaussian is flat at x = 0. 
Hence \\G,-N\\ X = 0{t- l l 2 ). 

This calculation shows that the fastest possible convergence 
is a power law with exponent a = 1/2, and provides an intu- 
itive reason why this is the case. If the rescaled shape function 
fit converges to a gaussian shape at a rate slower than f -1 / 2 , 
then the overall rate of convergence a could be slower than 
1/2. However, the numerical results in Sec. V B show that the 
rate is close to 1/2. We remark that for an observable which 
is not so intimately related to the geometrical structure of the 
lattice, the fine structure will in general be more complicated, 
and the above argument may no longer hold. 



VI. MAXWELLIAN VELOCITY DISTRIBUTION 

In this section we consider the effect of a non-constant dis- 
tribution of particle speeds [45]. A Maxwellian (gaussian) ve- 
locity distribution was used in polygonal and Lorentz channels 
in [37] and [20], respectively, in connection with heat con- 
duction studies. The mean square displacement was observed 
to grow asymptotically linearly, but the relationship with the 
unit speed situation was not discussed. A more complicated 
Lorentz gas with a gaussian distribution was studied in [38]. 

We show that the mean square displacement grows asymp- 
totically linearly in time with the same diffusion coefficient as 
for the unit speed case, but that the limiting position distribu- 
tion may be non-gaussian. For brevity we refer only to the 
position distribution throughout this section; the displacement 
distribution is similar. 



A. Mean square displacement 

Consider a particle located initially at (xo,vo), where vo 
has unit speed. Changing the speed of the particle does not 
change the path it follows, but only the distance along the path 
traveled in a given time. Denoting by 4> r v (xo,vo) the billiard 
flow with speed v starting from xo and with initial velocity in 
the direction of the unit vector vq, we have 



4>' v (xo,v ) =4> vr (x ,vo), 



(6.1) 



where the flow on the right hand side is the original unit-speed 
flow. If all speeds are equal to v, then the radially symmetric 
2D position probability density after a long time t is thus 



^^l V ) = 4^ eXP V ^D^, 



giving a radial density 

Pt ( r \ v ) = 



exp 



2Dvt \ADvt J 



(6.2) 



(6.3) 



(Throughout this calculation we neglect any fine structure.) 

If we now have a distribution of velocities with density 
pv(v), then the radial position density at distance r is 

f?\r)= r p t (r\v)p v (v)dv. (6.4) 

Jv=0 

The variance of the position distribution is then given by 

(x 2 )= f r 2 f^(r)dr (6.5) 

Jr=0 

r-ao 

= 4Dt v p v (v) dv =: 4Dtv, (6.6) 
Jo 

where v is the mean speed, after interchanging the integrals 
over r and v. 

We thus see that for any speed distribution having a finite 
mean, the variance of the position distribution, and hence the 
mean square displacement, grows asymptotically linearly with 
the same diffusion coefficient as for the uniform speed dis- 
tribution, having normalized such that v = 1 . We have veri- 
fied this numerically with a gaussian velocity distribution: the 
mean square displacement is indistinguishable from the unit 
speed case even after very short times [14]. 

B. Gaussian velocity distribution 

Henceforth attention is restricted to the case of a gaussian 
velocity distribution. For each initial condition, we gener- 
ate two independent normally-distributed random variables v\ 
and V2 with mean and variance 1 using the standard Box- 
Muller algorithm [18], and then multiply by a, which is a 
standard deviation calculated below. We use v\ and V2 as the 
components of the velocity vector v, whose probability den- 
sity is hence given by 



p{v) = p(v u v 2 ) = 



e -v2/2a 2 



-v 2 /2a 2 



a^/2% a^/2% ' 2na 2 



(6.7) 
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where v : 



speed v thus has density 



v 2 is the speed of the particle. The 



(6.8) 



and mean v = Oy/n/2. To compare with the unit speed dis- 
tribution we require v = 1, and hence a = \j2j%. As before, 
we distribute the initial positions uniformly with respect to 
Lebesgue measure in the billiard domain Q. 



C. Shape of limiting distribution 

The position density (6.4) is a function of time. However, 
the gaussian assumption used to derive that equation is valid 
in the limit when t — ► °°, so the central limit theorem rescaling 



fr d (r) := VF/, rad (^) 



(6.9) 



eliminates the time dependence in (6.4), giving the following 
shape for the limiting radial density: 



/ ra V) 



4D J V=Q 



exp 



ADv 



~4~ 



dv=: 



4D 



/, (6.10) 



denoting the integral by /. Mathematica [39] can evaluate 
this integral explicitly in terms of the Meijer G- function [40]: 

,.4 



I = G 



[ 256D 2 



0,0 



(6.11) 



See [41] and references therein for a review of the use of such 
special functions in anomalous diffusion. 

We can, however, obtain an asymptotic approximation to / 
from its definition as an integral, without using any properties 

J. 2 

of special functions, as follows. Define K(y) := tjj- + the 
negative of the argument of the exponential in (6.10). Then K 
has a unique minimum at v m j n := (r 2 / (2nD)) l l^ and we ex- 
pect the integral to be dominated by the neighborhood of this 
minimum. However, the use of standard asymptotic methods 



0, 



tends to 0, a 



is complicated by the fact that as r 
boundary of the integration domain. 

To overcome this, we change variables to fix the minimum 
away from the domain boundaries, setting w ;= v/v m ; n . Then 



/ = vv 



e -aL{ w ) . 



w=0 



(6.12) 



where a := 



and L(w) 



1 > 

1 with a minimum at 



2 V / vv 2 

Vfmin = 1- Laplace's method (see e.g. [42]) can now be ap 
plied, giving the asymptotic approximation 



(6.13) 



1 ^ Vinin c 



y/aL"(Wrmn) V3 



valid for large a, i.e. for large r. 
Hence 



r 4/3 



where 



C: 



2D^ H 2\32D 2 J 



(6.14) 



(6.15) 




FIG. 13: (Color online) The radial density function f t compared 
to the numerically calculated radial fine structure function rad , 
rescaled to converge to 1 and then vertically shifted for clarity. The 
demodulated radial density p, md is also shown, r = 2.3; b = 0.5; 
t = 100. 



D. Comparison with numerical results 

Figure 13 shows the numerical radial position density 
// ad (r) for a particular choice of geometrical parameters. We 
wish to demodulate this as in Sec. IV to extract the slowly- 
varying shape function, which we can then compare to the 
analytical calculation. 

The radial fine structure function md (r) must be calculated 
numerically, since no analytical expression is available. We do 
this by distributing 10 5 points uniformly on a circle of radius 
r and calculating the proportion of points not falling inside 
any scatterer. This we normalize so that rad (r) — > 1 as r — > 
°°, using the fact that when r is large, the density inside the 
circle of radius r converges to the ratio [r 2 — 7t(a 2 + b 2 )}/r 2 of 
available area per unit cell to total area per unit cell. We can 

u'ad 



then demodulate f, by <j) , setting 



fr(r) 

<j) md (ry/t)' 



(6.16) 



Figure 13 shows the demodulated radial density p r md (r) at 
two times compared to the exact solution (6. 10)— (6. 11), the 
asymptotic approximation (6.14)-(6.15), and the radial gaus- 
sian 275e _r2 / 2D . The asymptotic approximation agrees well 
with the exact solution except at the peak, while the numeri- 
cally determined demodulated densities agree with the exact 
long-time solution over the whole range of r. All three differ 
significantly from the gaussian, even in the tails. We conclude 
that the radial position distribution is non-gaussian. A similar 
calculation could be done for the radial displacement distribu- 
tion, but a numerical integration would be required to evaluate 
the relevant fine structure function. 

An explanation of the non-gaussian shape comes by con- 
sidering slow particles which remain close to the origin for 
a long time, and fast particles which can travel further than 
those with unit speed. The combined effect skews the re- 
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FIG. 14: (Color online) Comparison of the demodulated radial den- 
sity p, rad with the exact Meijer-G representation, the large-r asymp- 
totic approximation, and the radial gaussian with variance 2D. 



suiting distribution in a way which depends on the relative 
weights of slow and fast particles. 
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FIG. 15: (Color online) Rescaled ID marginal of the displacement 
density g and the demodulated version rj compared to the gaussian 
with variance 2D and to the asymptotic expression. The latter is not 
shown close to x = 0, where it drops to 0. r = 2.3; b = 0.5. 



VII. POLYGONAL BILLIARD CHANNEL 



E. ID marginal 

The ID marginal in the x-direction is shown in Fig. 15. 
Again there is a significant deviation of the demodulated den- 
sity from a gaussian. From (6.14), the 2D density at (x,y) is 
asymptotically 



C 



f(x,y) ~ ^ ex P 



-J3(^+/) 2 / 3 
from which the ID marginal f(x) is obtained by 



f(x,y)dy. 



(6.17) 



(6.18) 



It does not seem to be possible to perform this integration ex- 
plicitly for either the asymptotic expression (6.17) or the cor- 
responding exact solution in terms of the Meijer G-function. 
Instead we perform another asymptotic approximation start- 
ing, from the asymptotic expression (6.17). Changing vari- 
ables in (6.18) to z :=y/x and using the evenness in y gives 



where K : 



C\x\ 
2k 



exp 



2N2/3 



dz. 



(6.19) 



•|4/3 



Laplace's method then gives 



Cy/3 



IH 4/3 



(6.20) 



valid for large x. This is also shown in Fig. 15. Due to the 



|V3 



factor, the behavior near x = is wrong, but in the tails 



there is reasonably good agreement with the numerical results. 



In this section, we apply the previous ideas to a polygo- 
nal billiard channel. Polygonal models differ from Lorentz 
gases in that they are not chaotic in the standard sense, since 
the Kolmogorov-Sinai entropy and all Lyapunov exponents 
are zero due to the weak nature of the scattering from the 
polygonal sides [43]. Other indicators of the complexity of 
the dynamics of such systems are required: see e.g. [43] and 
references therein for a recent example. 

As far as we aware, there are few rigorous results on ergodic 
and statistical properties of these models [26, 44]. However, 
certain polygonal channels have been found numerically to 
show normal diffusion, in the sense that our property (a) is 
satisfied, i.e. the mean squared displacement grows asymptot- 
ically linearly: see e.g. [26, 37]. No convincing evidence has 
so far been given, however, that property (b), the central limit 
theorem, can be satisfied, although it was shown in [23] that 
(c) is satisfied for some random polygonal billiard models. 
Here we show that polygonal billiards can satisfy the central 
limit theorem. 



A. Polygonal billiard channel model 

We study a polygonal billiard introduced in [26]. The ge- 
ometry is shown in Fig. 16(a) and the channel in Fig. 16(b). 
We fix the angles (j>\ and 02 an d choose d such that the width 
of the bottom triangles is half that of the top triangle. This 
determines the ratio of h\ to I12 in terms of the angles (j>\ and 
02- We then require the inward-pointing vertices of each tri- 
angle to lie on the same horizontal line in order to prevent 
infinite horizon trajectories, giving h\-Yh% = h= \ and d = 
/;/(tan0i + j tan ^2), with /ij =dtan(j)\ and hi = (d/2) tan 02- 
We remark that in [26] it was stated that the area \Q\ = dh 
of the billiard domain is independent of <p2 when tj>\ is fixed, 
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(b) 



FIG. 16: (a) The geometry of the polygonal billiard unit cell, shown 
to scale with (fa = 7l/(2e). (b) Part of the polygonal channel with the 
same parameters. 



but this is not correct, since the expression for d shows that it 
depends on <fa, and we have fixed h = l. 

In [26] the parameters 0i = ^(^5— l)/8 and (fa = njq 
were used, with q € N and 3 < q < 9. For q > 5 normal 
diffusion was found, whereas for q = 3,4 it was found that 
(Ax 2 ) f ~ t a with a 7^ 1, so that property (a) is no longer satis- 
fied and we have anomalous diffusion. As far as we are aware, 
there is as yet no physical or geometrical explanation for this 
observed anomalous behavior, although presumably number- 
theoretic properties of the angles are relevant. 

We use the same (j)\ , but a value of 02 which is irrationally 
related to %, namely (fa = n/(2e) ~ 7C/5.44 (where e is the 
base of natural logarithms), since there is evidence that mix- 
ing properties are stronger for such irrational polygons [44]. 
In this case we find (Ax 2 ) t ~ t L008 , which we regard as asymp- 
totically linear, so that property (a) is again satisfied, with 
D = 0.3796 ±0.0009. 



B. Fine structure 

The shape of the displacement density was considered in 
[26] using histograms, but the results were not conclusive. 
Here we use our more refined methods to study the fine struc- 
ture of position and displacement distributions and to show 
their asymptotic normality. 

Figure 17 shows a representative position density ft(x). 
Following the method of Sec. IV D, we calculate the fine 
structure function h(x) as the normalized height of available 
space at position x; this is shown in the inset. We demodu- 
late ft by dividing by h to yield p,, which is again close to the 
gaussian with variance 2Dt . 

With the same notation as in Sec. IV F, we can also calcu- 
late the fine structure function of the displacement density. 
Taking the origin in the center of the unit cell in Fig. 16(a), we 




X 



FIG. 17: (Color online) Position density at t = 50 in the polygonal 
model with (fa = 7l/(2e). The inset shows h(x) over two periods. 



have 

h{x) = r^r (*tan0i + \x- ~d\ tan0 2 ) (7.1) 

for < x < d, with h being an even function and having pe- 
riod 2d. (The factor of 2d in (7.1) makes h a density per unit 
length.) The Fourier coefficients are h(Q) — 1 and 

6 <*> = ^/>K^) = i^ w a2) 

for ky^Q, where for m E Z we have 

!4tan(0i), if k odd 
8tan(0 2 ), iffc = 4m + 2 (7.3) 
0, if k = Am. 

C. Central limit theorem 

As for the Lorentz gas, we rescale the densities and distri- 
bution functions by \/t to study the convergence to a possible 
limiting distribution. Again we find oscillation on a finer and 
finer scale and weak convergence to a normal distribution: see 
Fig. 18. Figure 19 shows the time evolution of the demodu- 
lated densities f] t . There is an unexpected peak in the densities 
near x = for small times, indicating some kind of trapping 
effect; this appears to relax in the long time limit. Again we 
conjecture that we have uniform convergence of the demodu- 
lated densities fj to a gaussian density. 

Figure 20 shows the distance of the rescaled distribution 
functions from the limiting normal distribution, analogously 
to Fig. 12, for several values of (fa for which the mean square 
displacement is asymptotically linear. The straight line fitted 
to the graph for <p2 = 7C/(2e) has slope —0.212, so that the 
rate of convergence for this polygonal model is substantially 
slower than that for the Lorentz gas, presumably due to the 
slower rate of mixing in this system. A similar rate of decay 
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FIG. 18: (Color online) Rescaled displacement denstities compared 
to the gaussian with variance 2D. The inset shows the function <j> for 
this geometry. 




x 



FIG. 19: (Color online) Demodulated densities f\, for t = 100, t = 
1000 and t = 10000, compared to a gaussian with variance 2D. The 
inset shows a detailed view of the peak near x = 0. 

is found for fa = it /I, whilst fa = it/6 and fa = it/9 appear 
to have a slower decay rate. Nonetheless, the distance does 
appear to converge to for all these values of fa, providing 
evidence that the distributions are asymptotically normal, i.e. 
that the central limit theorem is satisfied. 

We remark that these convergence rate considerations will 
be affected if we have not reached the asymptotic regime, 
which would lead to an incorrect determination of the rele- 
vant limiting growth exponent and/or diffusion coefficient. 



VIII. CONCLUSIONS 

We have studied deterministic diffusion in diffusive bil- 
liards in terms of the central limit theorem. In a 2D periodic 
Lorentz gas model, where the central limit theorem is proved, 
we have shown that it is possible to understand analytically 
the fine structure occurring in the finite-time marginal posi- 
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FIG. 20: (Color online) Distance of the rescaled distribution func- 
tions from the limiting normal distribution for the polygonal model 
with different values of (j>2- The straight line is a fit to the large-time 
decay of the irrational case <j>2 = 7t/(2e). 



tion and displacement distribution functions, in terms of the 
geometry of a unit cell. Demodulating the observed densities 
by the fine structure allowed us to obtain information about 
the large-scale shape of the densities which would remain ob- 
scured without this demodulation: we showed that the demod- 
ulated densities are close to gaussian. 

We then studied the manner and rate of convergence to the 
limiting normal distribution required by the central limit the- 
orem. We were able to obtain a simple estimate of the rate of 
convergence in terms of the fine structure of the distribution 
functions. The demodulated densities appear to converge uni- 
formly to gaussian densities, which is a strengthening of the 
usual central limit theorem. 

We showed that imposing a Maxwellian velocity distribu- 
tion does not change the growth of the mean square displace- 
ment, but alters the shape of the limiting position distribution 
to a non-gaussian one. 

Finally we showed that similar methods can be applied to 
a polygonal billiard channel where few rigorous results are 
available, showing that the central limit theorem can be satis- 
fied by such models, but finding a slower rate of convergence 
than for the Lorentz gas. 

We believe that our analysis may have implications for the 
escape rate formalism for calculating transport coefficients 
(see e.g. [4]), where the diffusion equation with absorbing 
boundary conditions is used as a phenomenological model of 
the escape process from a finite length piece of a Lorentz gas: 
analyzing the fine structure in this situation could provide in- 
formation about the validity of this use of the diffusion equa- 
tion. We also intend to investigate models exhibiting anoma- 
lous diffusion using the methods presented in this paper. 
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